Stability of resonantly interacting heavy-light Fermi mixtures 
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We investigate a two- component mixture of resonantly interacting Fermi gases as a function of the 
ratio k of the heavy to the light mass of the two species. The diffusion Monte Carlo method is used 
to calculate the ground-state energy and the pair correlation function starting from two different 
guiding wave functions, which describe respectively the superfluid and the normal state of the gas. 
Results show that the mixture is stable and superfluid for mass ratios smaller than the critical value 
k c — 13 =b 1. For larger values of n simulations utilizing the wave function of the normal state are 
unstable towards cluster formation. The relevant cluster states driving the instability appear to be 
formed by one light particle and two or more heavy particles within distances on the order of the 
range of the interatomic potential. The small overlap between the wave function of the trimer bound 
state and the guiding wave function used to describe the superfluid state produces the unphysical 
stability of the superfluid gas above k c . 
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Two-component mixtures of fermionic gases with res- 
onant interspecies interactions constitute a fundamen- 
tal paradigm in the physics of ultracold gases. In the 
last decade a large amount of experimental and the- 
oretical work has been devoted to the study of mix- 
tures of two different hyperfine states in a single species 
gas (either 6 Li or 40 K) in the close vicinity of a Fano- 
Feshbach resonance [l| . At low enough temperature and 
for balanced populations of the two components, the 
gas has been found to be superfluid irrespective of the 
strength of interactions. The pairing mechanism, how- 
ever, changes as a function of coupling from the Cooper 
pairing in the Bardeen-Cooper-SchriefTer (BCS) region to 
the molecule formation in the Bose-Einstein condensate 
(BEC) regime. The unitary limit, corresponding to the 
resonance position, has attracted the greatest attention 
because of its universal character, as the inverse Fermi 
momentum kp 1 and the Fermi energy ep are the only 
relevant length and energy scale, respectively. 

More recent important experimental advances concern 
the realization of resonantly interacting mixtures with 
two different fermionic species 0-0 • In this case a new 
degree of freedom, namely the ratio k = M/m of the 
heavy to the light mass of the two species, enters the 
problem leading to interesting theoretical predictions for 
the superfluid ground state (9l4ll|. 

The effects of mass asymmetry between the heavy (h) 
and light (I) components of the mixture have been thor- 
oughly investigated at unitarity for systems of few atoms, 
both in free space and in harmonic traps EMI- It has 
been found [HI, [l3[ that an infinite number of Efimov 
states composed by two h and one / particles appear when 
the mass ratio n > 13.6. In the range 8.62 < n < 13.6 
interactions can be fine tuned to produce triplet bound 
states [14]]. Furthermore, few-body calculations of sys- 
tems composed by four and five particles indicate that 



a bound state can appear at n c± 10.4 for the system 
(3/i+lZ) and at k ~ 9.8 for (4A+ 1/) [13, 

Many-body heavy-light Fermi mixtures at unitarity 
have also been studied using quantum Monte Carlo 
(QMC) methods by von Stecher et al. [19] in harmonic 
traps and by Gezerlis et al [20| in homogeneous systems. 
The aim of Ref. [20| was mainly to analyze the particular 
case of 40 K- 6 Li mixtures with k ~ 6.5 as a function of 
the population imbalance. Interestingly, they comment 
that the normal state of the gas collapses to a many- 
body bound state (cluster) for mass ratios n ~ 12. This 
finding points to the existence at unitarity of a critical 
value k c above which the homogeneous gas is mechani- 
cally unstable against density fluctuations. The strongly 
interacting 40 K- 6 Li mixture, already produced in experi- 
ments |2r4Sl] , has a mass ratio too small for observing this 
instability and one can think that other possibilities, as 
for instance 173 Yb- 6 Li [2l|, could be brought to the res- 
onant regime in a near future. Moreover, even for lighter 
species one could also imagine that the mass ratio can 
be tuned in a nearly continuous way through proper con- 
finement in an optical lattice (in this case, one relies on 
effective masses) [10j . 

In the present work we investigate the stability of 
heavy-light Fermi mixtures at unitarity by carrying out 
a fully microscopic study using the diffusion Monte Carlo 
(DMC) method within the fixed-node (FN) approxima- 
tion. This technique has become nowadays a standard 
tool used extensively in many problems in condensed- 
matter physics including ultracold gases 22]. We use 
an approach similar to Ref. [23| to calculate the ground- 
state energy and the pair correlation function starting 
from a guiding wave function that describes either the 
superfluid or the normal state of the gas. We find that 
the latter shows an instability towards cluster formation 
at the critical mass ratio k c = 13 ± 1. In the case of 
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the superfluid state, the apparent lack of instability can 
be understood as an artifact in the nodal surface of the 
corresponding guiding wave function. 

The Hamiltonian of a mixture of heavy and light 
fermions with masses M and m is written as 



fc2 N h h 2 
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where denotes the Laplacian operator for the heavy 

(light) particles. Periodic boundary conditions in a cu- 
bic box of volume V = L 3 are used and we restrict our 
study to balanced populations, Nh = Ni = N/2, N being 
the total number of particles. Intraspecies interactions 
are neglected while short-range interspecies interactions 
are modeled by a square- well potential with depth Vq: 
V(r) = — Vo if r < R and V(r) = elsewhere. For this 
potential the s-wave scattering length is known analyti- 
cally: a = R [1 - tfm(K R)/(K R)], with K 2 = fiV /H 2 
and \i = (M _1 -f m -1 ) -1 the reduced mass. Our interest 
is in the unitary limit where \a\ = oo and the resonance 
condition gives in our case KqR = 7r/2, corresponding 
to the appearance of the first two-body bound state in 
the well. To ensure universality the value of R is chosen 
much smaller than the mean interparticle distance deter- 
mined by the density n = N/V: as in Ref. [23| we use 
nR 3 = KT 6 . 

We carry out simulations using two different guiding 
wave functions. The first is a Jastrow-Slater (JS) func- 
tion given by 
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where f(r) is a non- negative function of the interparti- 
cle distance reproducing at short separations the zero- 
energy solution of the two-body problem with the poten- 
tial V(r) and satisfying the boundary condition f'(r = 
L/2) = for its derivative at half the box width. The 
pair orbital in the determinant is given by the combi- 
nation of plane waves ip(r) = Xln 2 <4 cos (k n * r), where 
k n = 2it(n Xl n yi n z ) / L are the wave vectors determined 
by the integers n XjVjZ = 0, ±1, ... and the sum is restricted 



over the filled shells n 
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< 4, correspond- 



ing to TV = 66 particles used in the simulations. The 
determinant in (|2j) fixes the nodal surface of the many- 
body wave function to that of a non-interacting gas, be- 
ing therefore incompatible with off-diagonal long-range 
order. As a consequence, the Jastrow-Slater wave func- 
tion (|2|) describes the normal state of the gas. Simula- 
tions of the superfluid state are instead performed using 
the superfluid wave function 



* T (n,...,rjv) =det[0(|rf -rj 
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where the pair orbital <j)(r) is chosen as the two-body so- 
lution, similarly to the Jastrow factor f(r) in Eq. (j2]). 




FIG. 1. (color online). Energy per particle in units of cifg 
as a function of the mass ratio k. Diamonds correspond to 
the results for the superfluid state. Circles and squares cor- 
respond instead to simulations of the normal state with a 
number of walkers N w = 1000 and 5000, respectively. The 
insets show for this latter case a two-dimensional projected 
density for two different values of k. 



In fact, Eq. (j3j) corresponds to a BCS wave function 
characterized by a spherically symmetric order param- 
eter. In the FN-DMC algorithm an important system- 
atic bias concerns the number of walkers N w , i. e. points 
R = (n, rjv) in configuration space where the wave 
function is sampled. The FN-DMC energies should be 
estimated in the limit N w — » oo. However, if the guid- 
ing function is reasonably accurate, the dependence 
on N w can be neglected for some large enough number 
of walkers. On the contrary, a strong dependence on 
N w points out the effect of some relevant physics not 
accounted for by ^t- 

Mean-field theory has been applied to the problem of 
the BCS-BEC crossover at T = with mass asymme- 
try [9| , but the results obtained show that both the chem- 
ical potential and the gap parameter, once the Fermi en- 
ergy is rescaled with the reduced mass ep = h 2 k 2 F /An 
where kp = (37r 2 n) 2 / 3 , are unchanged compared to 
the corresponding quantities of the equal mass problem. 
Non-trivial effects arising from the mass ratio appear in 
the expression of the gap parameter if one includes per- 
turbative corrections following the approach by Gorkov 
and Melik-Barkhudarov [10|. 

We calculate the energy as a function of the mass 
ratio up to large values of k, in particular larger than 
n = 13.6 where bound trimers are expected to appear. 
The FN-DMC energies of the superfluid and normal state 
are shown in Fig. [TJ For n — 1 we find that the su- 
perfluid state is lower in energy than the normal state: 
i = 0.42 ± 0.01 vs. i = 0.54 ± 0.01 • Here f is the 

proportionality coefficient between the energy per parti- 
cle of the interacting system E/N = ^cifg and the one 
of a free Fermi gas cifg = 3ep/5. For n > 1 both en- 
ergies decrease slightly with increasing in agreement 
with the findings of Ref. [20] reporting for k = 6.5 a 
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FIG. 2. (color online): Energy per particle of the normal 
state with N = 14 as a function of the inverse of the number 
of walkers N w for k = 10, 12, and 14 (respectively red, green 
and blue symbols). The inset shows the energy probability 
distribution for the same values of k and N w = 800. 
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FIG. 3. (color online): Pure estimators for the two-body 
radial distribution functions #2 h (r) and g^i?) for two different 
values of k. Distances are in units of the interaction range R 
and in these units k^ 1 ~ 32. The pair correlation function of 
an ideal Fermi gas is also shown for comparison. 



~ 5% decrease of the energy of the superfluid state com- 
pared to the equal mass case. We notice that in the range 
ft < 13 the energy of the normal state lies above the one 
of the superfluid state implying that the latter is the sta- 
ble ground-state of the system. For ft > 14 the behavior 
in the two cases is completely different: the energy calcu- 
lated using the superfluid guiding wave function slightly 
decreases with ft, whereas the energy obtained using the 
normal gas wave function shows a fast-growing instability 
towards negative energies. According to the variational 
nature of fixed- node calculations, we conclude that the 
gas becomes unstable against cluster formation around 
ft ~ 13. This instability is explicitly observed in the par- 
ticle configurations of the FN-DMC simulation. In the 
insets of Fig. [TJ we show typical snapshots of the parti- 
cle positions for two mass ratios, below and above the 
instability. For ft = 10 the gas is structureless and ho- 
mogeneous, whereas for ft = 14 the system has formed 
clusters indicating a spinodal decomposition. 

The dependence of the energy of the normal state on 
the number of walkers is first shown as a function of the 
mass ratio in Fig. [J for N w = 1000 and 5000 and is fur- 
ther analyzed in Fig. [2j For ft = 10 and 12 the energy 
shows a linear decrease with increasing N w towards the 
limit 1/N W = with a well defined asymptotic value. On 
the contrary, for ft = 14 the energy decreases fast with 
increasing N w and for large N w we observe negative ener- 
gies corresponding to the formation of self-bound many- 
body clusters. The appearance of clusters is also reflected 
in the large increase of the statistical error when ft = 14 
and N w is large. Additional insight on the increase in 
variance can be inferred from the probability distribu- 
tion of local energies. The energy histograms are shown 
in the inset of Fig. O For ft = 10 and 12 there is a 
well defined Gaussian peak, centered on the mean value, 
whose width determines the variance of the energy esti- 



mate. For ft = 14 the picture is appreciably different: 
large tails in the distribution are developed, mainly to- 
wards lower energies. 

The onset of instability can also be observed in the 
two-body radial distribution function of the gas. The 
functions #2^( r )> with a, f3 = h, I are proportional to 
the probability of finding two particles of a and j3 type 
at distance r. Their behavior at unitarity in the case 
ft = 1 has been analyzed in Refs. [23|, \2j\. For like par- 
ticles #2 a ( r ) i s very similar to that of an ideal Fermi 
gas, because Pauli exclusion principle largely overcomes 
the induced correlations mediated by interactions with 
the other component. Instead a large peak is formed in 
(r) at short distances for unlike particles due to the 
their strong attractive interactions. A similar behavior is 
expected in ^^{r) when ft > 1. The tendency of cluster 
formation is instead visible in the distribution function 
of like particles. In Fig. [3j we report results for #2 h ( r ) 
and #2 (r) obtained for mass ratios ft = 10 and 14. The 
distribution function of light particles is insensitive to 
the value of ft remaining close to the ideal gas result. 
On the contrary, #2 h ( r ) shows a large peak at short dis- 
tances for ft = 14, indicating the formation of bound 
states involving one light particle. The Pauli exclusion 
principle brings eventually g^ h down to zero on length 
scales smaller than the interaction range R. 

An estimate of the critical mass ratio ft c characterizing 
the onset of the cluster instability can be obtained by 
integrating the radial distribution function up to a cut 
off length R C1 
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for the hh and 11 distribution functions. The value of 
R c is chosen to be large compared to the range R of the 
interatomic potential, giving the typical size of cluster 
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FIG. 4. (color online): Results for the integrated distribution 
function g 1 2 lt (R c ) denned in Eq. (g) as a function of the mass 
ratio k for the case of 11 and hh particles. In both cases, 
R c — SR. 




FIG. 5. Probability distribution as a function of the distances 
| ri — r l \ and \r% — r l \ for a three-particle system composed 
by two heavy and one light particles with mass ratio k = 20. 
(a): normal wave function, (b): superfluid wave function. 



r l \) sin(27r(xj' — x^j/L), corresponding to the JS normal 
wave function (O, and ^ T = </>(l r i ~ r 1) ~ <Kl r 2 ~ r 1)> 
corresponding to the superfluid wave function (|3]). The 
nodal surfaces bare important differences: the molecu- 
lar orbital of the superfluid wave function is spherically 
symmetric while the one of the JS function possesses 
the cubic symmetry of plane waves in a box. Forma- 
tion of a trimer bound state is strongly suppressed in the 
case of the superfluid wave function. In fact, the two 
heavy particles repel each other due to the Pauli prin- 
ciple and prefer to stay apart, while the light particle 
attracts both heavy particles. If one assumes that in 
the trimer the heavy particles will preferentially stay at 
the same distance from the light particle, then a node 
appears in when \r± — r l \ = \r^ — r l \. On the con- 
trary, in the JS guiding function there is an enhanced 
probability of arranging two heavy particles at the same 
distance from the light particle. This different behav- 
ior of the two wave functions is clearly shown in Fig. [5] 
where we show the spatial dependence of the probabil- 
ity of finding the two heavy particles at the distances 
\r± — r l \ and \r^ — r l \ from the light particle. This cor- 
relation function vanishes when both distances approach 
zero according to the Pauli exclusion principle acting on 
the heavy particles. At the mass ratio k = 20 shown in 
the figure, one sees that the probability of finding the 
three particles close to each other is greatly increased in 
the case of the JS function. The peak corresponds to 
the sum \r^ — r l \ + \r% — r l \ ~ R. On the contrary, for 
the superfluid wave function the probability of finding 
— r l \ = | rJj — r l \ is greatly suppressed and the overlap 
of this wave function with cluster states composed of two 
heavy and one light particle is therefore very small. 



states, but small compared to the mean interparticle dis- 
tance. In Fig. |4]we show the integrated values of #2 h ( r ) 
and g l 2(r), with R c = 8R, as a function of the mass ra- 
tio. For the 11 distribution function the results of g 2 nt are 
essentially independent of ft, in agreement with the be- 
havior shown in Fig. [3j On the contrary, in the hh case, 
the corresponding g 2 nt increases linearly with the mass 
ratio up to k ~ 11. For k > 13 the increase is still linear, 
but the slope is significantly larger. In Fig.|4]we show the 
linear fit to the data in the regime k < 11 and k > 13. 
We checked that by changing the value of R Cl provided 
it remains within the range discussed above, the quali- 
tative picture shown in Fig. |4] remains the same and in 
particular the value of n corresponding to the change of 
slope in g2 lt (R c ) for the hh particles. 

An important question to understand is the apparent 
stability of simulations carried out using the superfluid 
guiding function (J3j) when the mass ratio is larger than n c . 
In order to gain insight on that, we consider a simple sys- 
tem composed by two heavy and one light particles. We 
perform FN-DMC simulations of this three-body prob- 
lem using as guiding functions = Y\i=i 2 f(\ T i ~ 



In conclusion, we investigate the stability at unitarity 
of a many-body mixture of two fermionic species with 
mass ratio n > 1 using the FN-DMC method. Our results 
show that the mixture is superfluid and stable up to a 
critical mass ratio n c — 13 ± 1. For larger mass ratios, an 
instability sets in towards the formation of clusters. The 
study of the pair correlation function indicates that the 
relevant cluster states driving the instability are formed 
by one light particle and two or more heavy particles 
within distances on the order of the range R of the po- 
tential. Within our approach, no many-body instability 
is observed for k < 12 where four and five-body bound 
states have been found in few-body calculations [13, [3 • 
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